\documentclass[10pt]{article}
\usepackage{amsmath}
\usepackage{graphicx}

\newcommand{\problem}[1]{\noindent{\sc Problem #1.}}
\newcommand{\solution}{\noindent{\sc Solution.}}
\newcommand{\complexity}{\noindent{\sc Complexity. }}
\newcommand{\answer}{\noindent{\sc Answer. }}

\setlength{\parskip}{1ex plus 0.5ex minus 0.2ex}

\begin{document}

\problem{339}
Peredur vab Efrawc.

\emph{
``And he came towards a valley, through which ran a river; and the borders of the valley were wooded, and on each side of the river were level meadows. And on one side of the river he saw a flock of white sheep, and on the other a flock of black sheep. And whenever one of the white sheep bleated, one of the black sheep would cross over and become white; and when one of the black sheep bleated, one of the white sheep would cross over and become black."} -- Peredur the Son of Evrawc

Initially each flock consists of $n$ sheep. Each sheep (regardless of color) is equally likely to be the next sheep to bleat. After a sheep has bleated and a sheep from the other flock has crossed over, Peredur may remove a number of white sheep in order to maximize the expected final number of black sheep. Let $E(n)$ be the expected final number of black sheep if Peredur uses an optimal strategy.

You are given that $E(5) = 6.871346$ rounded to 6 places behind the decimal point.

Find $E(10000)$ and give your answer rounded to 6 places behind the decimal point.

\solution

Starting with a configuration of $b$ black sheep and $w$ white sheep, let $f(b,w)$ be the expected number of black sheep in the end where Peredur has the chance to remove a certain number of white sheep before any sheep bleats, and $e(b,w)$ be the expected number of black sheep in the end where some sheep bleats first before Peredur does anything. It is easy to find that
\begin{align}
f(b,w) &= \max_{0 \le k \le w} e(b, k) \notag \\
&= \max \left\{ e(b,w), \, \max_{0 \le k \le w-1} e(b, k) \right\} \notag \\
&= \max \{ e(b,w), f(b,w-1) \}
\end{align}
We also note that
\[
e(b,w) = \frac{b}{b+w} f(b+1,w-1) + \frac{w}{b+w} f(b-1, w+1)
\]
for $b>0, w>0$, and $e(b,0)=b, e(0,w)=0$. Combining these equations, we get
\[
f(b, w) = \max \left\{ f(b,w-1), \, \frac{b}{b+w} f(b+1,w-1) + \frac{w}{b+w} f(b-1, w+1) \right\}
\]
for $b>0, w>0$, and $f(b,0) = b, f(0,w) = 0$.

%Start with the simple cases. When there is one sheep, we have
%\begin{align}
%f(1,0) &= 1 \notag \\
%f(0,1) &= 0 \notag
%\end{align}
%When there are two sheep,
%\begin{align}
%f(2,0) &= 2 \notag \\
%f(1,1) &= \max \left\{ f(1,0), \frac12 f(2,0)+\frac12 f(0,2) \right\} = 1 \notag \\
%f(0,2) &= 0 \notag
%\end{align}
%When there are three sheep,
%\begin{align}
%f(0,3) &= 0 \notag \\
%f(1,2) &= \max \left\{ f(1,1), \frac13 f(2,1) + \frac23 f(0,3) \right\} = 1 \notag \\
%f(2,1) &= \max \left\{ f(2,0), \frac23 f(3,0) + \frac13 f(1,2) \right\} = \frac73 \notag \\
%f(3,0) &= 3 \notag
%\end{align}
%When there are four sheep,
%\[
%f(4,0)=4, f(3,1)=\frac{43}{12}, f(2,2)=\frac{7}{3}, f(1,3)=1, f(0,4)=0
%\]
%When there are five sheep,
%\[
%f(5,0)=5, f(4,1)=\frac{157}{33}, f(3,2)=\frac{125}{33}, f(2,3)=\frac{7}{3},
%f(1,4)=1, f(0,5)=0
%\]

%
%In general, let's revisit equation
%\[
%f(b, w) = \max \left\{ f(b,w-1), \, \frac{b}{b+w} f(b+1,w-1) + \frac{w}{b+w} f(b-1, w+1) \right\}
%\]

To solve this system of equations, let $n = b + w$ be the total number of sheep. For a given $n$, let $f_n(k) \equiv f(k, n - k)$, and we can rewrite the equation as
\[
f_n(k) = \max \left\{ f_{n-1}(k), \, \frac{k}{n} \times f_n(k+1) + \left(1 - \frac{k}{n}\right) \times f_n(k-1) \right\}
\]
and $f_n(0) = 0, f_n(n) = n$. Below is a plot of the values of $f_n(k)$ for $1 \le n \le 10$.
\begin{center}
\includegraphics[height=45mm]{p339.png}
\end{center}

We observe (without proof) that
\[
f_n(k) = \left\{
\begin{array}{ll}
f_{n-1}(k) & \text{for } 1 \le k \le \lfloor n/2 \rfloor, \\
\frac{k}{n} \times f_n(k+1) + \left(1 - \frac{k}{n}\right) \times f_n(k-1) & \text{for } \lfloor n/2 \rfloor < k < n. \\
\end{array} \right.
\]
Once we have the value of $f_{n-1}(\lfloor n/2 \rfloor)$, we can solve the linear equation system of $f_n(k)$ for $\lfloor n/2 \rfloor < k < n$.

To solve the linear system, we use Gaussian elimination. Let
\[
f_n(k) = p_k f_n(k+1) + q_k
\]
where $p_k$ and $q_k$ are coefficients to be determined. By comparing with the previous equation, we find
\[
p_k = \frac{1-a_k}{1-a_k p_{k-1}}, \, q_k = \frac{a_k q_{k-1}}{1-a_k p_{k-1}}
\]
for $i < k < j$, with starting value $p_i = 0, q_i = x'_i$. Since we already know $x_j = x'_j$, we can find all the $x_k$'s backwards.


%%To obtain a solution to the above equations, we instead solve the following system of inequalities
%%\[
%%x_k \ge b_k, x_k \ge a_k x_{k-1} + (1-a_k) x_{k+1}
%%\]
%%to minimize each $x_k$. To proceed, we use an approach similar to Gaussian elimination. Suppose we have the following relation
%%\[
%%x_k \ge p_k x_{k+1} + q_k
%%\]
%%where $p_k$ and $q_k$ are coefficients to be determined. Because $x_0 = 0$, we have $p_0 = 0, q_0 = 0$. For $1 \le k \le n-1$, substitute
%%\[
%%x_{k-1} \ge p_{k-1} x_k + q_{k-1}
%%\]
%%into the above equation, and we get
%%\[
%%x_k \ge a_k ( p_{k-1} x_k + q_{k-1} ) + (1-a_k) x_{k+1}
%%\]
%%Rearranging terms, we get
%%\[
%%(1-a_k p_{k-1}) x_k \ge a_k q_{k-1} + (1-a_k) x_{k+1}
%%\]
%%Since the right-hand side is greater than zero, and $x_k$ is greater than zero, so $1-a_k p_{k-1}>0$. Therefore
%%\[
%%x_k \ge \frac{1-a_k}{{1-a_k p_{k-1}}}  x_{k+1} + \frac{a_k q_{k-1}}{1-a_k p_{k-1}}
%%\]
%%This gives
%%\[
%%p_k = \frac{1-a_k}{{1-a_k p_{k-1}}}, \, q_k = \frac{a_k q_{k-1}}{1-a_k p_{k-1}}
%%\]
%%Since we know $x_n = n$, we can use the above equations to find out all the $x_k$ backwards from $x_{n-1}$ to $x_1$.
%%
%%We summarize the algorithm to find $E(N) = e(N,N)$ as follows.
%%
%%Step 1. Let $n = 1$.
%%
%%Step 2. For a given $n$, let $p_0=0, q_0=0$. Compute $p_k, q_k$ using the iterative formula
%%\[
%%p_k = \frac{1-a_k}{{1-a_k p_{k-1}}}, \, q_k = \frac{a_k q_{k-1}}{1-a_k p_{k-1}}
%%\]
%%for $1 \le k \le n-1$, where $a_k = (n-k)/n$.

%We use the following iterative approach to solve this system of linear equations with lower bound.
%
%Step 1. Let $x_k = b_k$ for every $1 \le k \le n-1$ as the initial solution and mark each $x_k$ \emph{not bumped}. In the following steps, $x_k$ will never decrease, so we needn't worry about the condition $x_k \ge b_k$ any more.
%
%Step 2. Given a solution $(x_1, \ldots, x_{n-1})$, check whether the equation
%\[
%x_k \ge a_k x_{k-1} + (1-a_k) x_{k+1}
%\]
%holds for each $k$ starting from $n-1$ backwards to $1$ where $x_k = b_k$ (i.e where $x_k$ is not bumped). If $x_k$ does not satisfy the condition, then mark $x_k$ as bumped, and set $x_k = a_k x_{k-1} + (1-a_k) x_{k+1}$. If all $x_k$'s satisfy the condition, then we have found the solution.
%
%Step 3. For each range $i < j$ where $x_i$ and $x_j$ are not bumped but all $x_k$ between them are bumped, (note that there can be multiple such ranges), solve the system of linear equations between $i$ and $j$:
%\[
%x_k = a_k x_{k-1} + (1-a_k) x_{k+1}
%\]
%where $i < k < j$, given $x_i = x'_i$ and $x_j = x'_j$. This system is easy to solve. Let
%\[
%x_k = p_k x_{k+1} + q_k
%\]
%where $p_k$ and $q_k$ are coefficients to be determined. By comparing with the previous equation, we find
%\[
%p_k = \frac{1-a_k}{1-a_k p_{k-1}}, \, q_k = \frac{a_k q_{k-1}}{1-a_k p_{k-1}}
%\]
%for $i < k < j$, with starting value $p_i = 0, q_i = x'_i$. Since we already know $x_j = x'_j$, we can find all the $x_k$'s backwards.
%
%Step 4. Go to step 2.
%
%This algorithm is guaranteed to finish because in each iteration, more variables become \emph{bumped}. The time complexity is $\mathcal{O}(n^2)$.

%%Note: We show that the system of equations
%%\[
%%x_k = \max \left\{ b_k, \, a_k x_{k-1} + (1-a_k) x_{k+1} \right\}
%%\]
%%is equivalent to the linear programming problem
%%\[
%%\min \sum_{k=1}^{n-1} x_k
%%\]
%%subject to
%%\[
%%x_k \ge b_k, \, x_k \ge a_k x_{k-1} + (1-a_k) x_{k+1},
%%\]
%%which can be solved by standard methods.
%%
%%Let $(x_1, \ldots, x_{n-1})$ be a solution to the linear programming problem. Then we must have $x_k = \max \{ b_k, \, a_k x_{k-1} + (1-a_k) x_{k+1} \}$. Suppose that this is not true, then there exists an $x_k$ such that
%%\[
%%x_k = \delta + \max \{ b_k, \, a_k x_{k-1} + (1-a_k) x_{k+1} \}
%%\]
%%with $\delta > 0$. Now replace $x_k$ with $(x_k - \delta)$ in the solution. We can see that all constraints are still satisfied, but the objective function is smaller. This contradicts with the precondition that $\{x_k\}$ is an optimal solution to the LP problem.




\answer
$e(10000,10000) = 19823.542204$

\complexity

Time complexity: $\mathcal{O}(n^2)$

Space complexity: $\mathcal{O}(n)$

\end{document} 